submitted to ApJ 



TIME-DEPENDENT STOCHASTIC PARTICLE 

ACCELERATION IN ASTROPHYSICAL PLASMAS: 

EXACT SOLUTIONS INCLUDING 

^ : MOMENTUM-DEPENDENT ESCAPE 

O 

^ : Peter A. Becker^'^ 

;h 

^ . Center for Earth Observing and Space Research, 

^ , George Mason University 

^ : Fairfax, VA 22030-4444, USA 



>- ■ Truong Le^ & Charles D. Dermer'* 

E. 0. Hulburt Center for Space Research 



O 



■^ ' Naval Research Laboratory, 

§ ■ Washington, DC 20375, USA 

o' 

Ph: abstract 
o 

c/3 . Stochastic acceleration of charged particles due to interactions with magne- 
tohydro dynamic (MHD) plasma waves is the dominant process leading to the 



1 



formation of the high-energy electron and ion distributions in a variety of astro- 
physical systems. Collisions with the waves influence both the energization and 
the spatial transport of the particles, and therefore it is important to treat these 
two aspects of the problem in a self-consistent manner. We solve the representa- 
tive Fokker-Planck equation to obtain a new, closed-form solution for the time- 
dependent Green's function describing the acceleration and escape of relativistic 
ions interacting with Alfven or fast-mode waves characterized by momentum dif- 
fusion coefficient D{p) oc p'^ and mean particle escape timescale tesc(p) oc p'^~^. 
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where p is the particle momentum and q is the power-law index of the MHD wave 
spectrum. In particular, we obtain solutions for the momentum distribution of 
the ions in the plasma and also for the momentum distribution of the escaping 
particles, which may form an energetic outflow. The general features of the so- 
lutions are illustrated via examples based on either a Kolmogorov or Kraichnan 
wave spectrum. The new expressions complement the results obtained by Park 
and Petrosian, who presented exact solutions for the hard-sphere scattering case 
(g = 2) in addition to other scenarios in which the escape timescale has a power- 
law dependence on the momentum. Our results have direct relevance for models 
of high-energy radiation and cosmic-ray production in astrophysical environments 
such as 7-ray bursts, active galaxies, and magnetized coronae around black holes. 
In particular, we outline an application of the new results to black-hole sources 
that produce outflows of relativistic hadrons, with associated predictions that 
can be tested using GLAST. 

Subject headings: acceleration of particles — methods: analytical — cosmic rays 
— black hole physics — plasmas — galaxies: jets 



1. INTRODUCTION 

Observations of high-energy radiation from a variety of astrophysical sources imply 
the presence of significant populations of nonthermal (often relativistic) particles. Much 
of the observational data is characterized by strong variability on very short timescales. 
Nonthermal distributions are naturally produced via the Fermi process, in which particles 
interact with scattering centers moving systematically and/or randomly. The first-order 
Fermi mechanism treats particle acceleration in converging fiows, such as shocks, and is 
thought to be important at the Earth's bow shock, in certain classes of solar fiares, for 
cosmic-ray acceleration by supernova remnants, and in sources with relativistic outfiows, 
including blazars and 7-ray bursts (for reviews, see Blandford & Eichler 1987; Kirk 1994). 
The second-order process, as it was originally conceived by Fermi, involved the stochastic 
acceleration of particles scattering with randomly moving magnetized clouds (Fermi 1949). 
Later refinements of this idea replaced the magnetized clouds with magnetohydrodynamic 
(MHD) waves (e.g., Melrose 1980). The second-order, stochastic Fermi process now finds 
application in a wide range of astrophysical settings, including solar fiares (Miller et al. 1990; 
Liu et al. 2004a), clusters of galaxies (Schlickeiser et al. 1987; Petrosian 2001; Brunetti et 
al. 2004), the Galactic center (Liu et al. 2004b; Atoyan & Dermer 2004), and 7-ray bursts 
(Waxman 1995; Dermer & Humi 2001). 
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The standard approach to modeling the acceleration of nonthermal particles via inter- 
actions with MHD waves involves obtaining the solution to a steady-state transport equa- 
tion that incorporates treatments of systematic and stochastic Fermi acceleration, radiative 
losses, and particle escape (e.g., Schlickeiser 1984; Schlickeiser & Steinacker 1989; Liu et al. 
2006). However, the prevalence of variability on short timescales in many sources calls into 
question the validity of the steady-state interpretation. Park & Petrosian (1995) provided 
a comprehensive review of the various time- dependent solutions that have been derived in 
the past 30 years. In most cases, it is assumed that the momentum diffusion coefficient, 
D[p), and the mean particle escape timescale, tesc{p), have power-law dependences on the 
particle momentum p. Although the set of existing solutions covers a broad range of possible 
values for the associated power-law indices, there are certain physically interesting cases for 
which no analytical solution is currently available. For example, Dermer et al. (1996) have 
shown that for the stochastic acceleration of relativistic particles due to resonant interactions 
with plasma waves in black-hole magnetospheres, one obtains D{p) oc p'^ and tcsc{p) oc p'^~'^, 
where q is the index of the wavenumber distribution (see § 2). The analytical solution for the 
time-dependent Green's function in this situation is of particular physical interest, but it has 
not appeared in the previous literature. This has motivated us to reexamine the associated 
Fokker-Planck transport equation for this case, and to obtain a new family of closed-form 
solutions for the secular Green's function. The resulting expression, describing the evolution 
of an initially monoenergetic particle distribution, complements the set of solutions discussed 
by Park & Petrosian (1995). Our primary goal in this paper is to present a detailed deriva- 
tion of the exact analytical solution and to demonstrate some of its key properties. Detailed 
applications of our results to the modeling of particle acceleration in black-hole accretion 
coronae and other astrophysical environments will be presented in a separate paper. 

The remainder of the paper is organized as follows. In § 2 we review the fundamental 
equations governing the acceleration of charged particles interacting with plasma waves. 
The transport equation is solved in § 3 to obtain the time-dependent Green's function, and 
illustrative results are presented in § 4. The astrophysical implications of our work are 
summarized in § 5, and additional mathematical details are provided in the Appendix. 



2. FUNDAMENTAL EQUATIONS 

Charged particles in turbulent astrophysical plasmas are expected to be accelerated via 
interactions with whistler, fast-mode, and Alfven waves propagating in the magnetized gas. 
Here we consider a simplified isotropic description of the wave energy distribution, denoted by 
W{k), where W{k)dk represents the energy density due to waves with wavenumber between 
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k and k + dk. The transport formalism we consider assumes a power-law distribution for the 
wave-turbulence spectrum, which implies definite relations between the momentum diffusion 
coefficient and the momentum-dependent escape timescale (e.g., Melrose 1974; Schlickeiser 
1989a,b; Dermer et al. 1996). MHD waves injected over a narrow range of wavenumber 
cascade to larger wavenumbers, forming a Kolmogorov or Kraichnan power spectrum over 
the inertial range with W{k) oc k~'^, where q = 5/3 and q = 3/2 for the Kolmogorov 
and Kraichnan cases, respectively (e.g., Zhou & Matthaeus 1990). The specific forms we will 
adopt for the transport coefficients in § 2.2 are based on the physics of the resonant scattering 
processes governing the wave-particle interactions (see Dermer et al. 1996). We assume that 
the nonthermal particle distribution is isotropic, and focus on a detailed treatment of the 
propagation of particles in momentum space due to wave-particle interactions. The spatial 
aspects of the transport (i.e., the confinement of the particles in the acceleration region) will 
be treated in an approximate manner using a momentum-dependent escape term. 



2.1. Transport Equation 

The fundamental transport equation describing the propagation of particles in the mo- 
mentum space can be written in the fiux-conservation form (e.g., Becker 1992; Schlickeiser 
1989a,b) 

— = — r— <!» A{p)f-D{p) — 



dt p^ dp \ 



' +^. (i; 



tcscip) 47rp2 

where p is the particle momentum, f{p,t) is the particle distribution function, D{p) denotes 
the momentum diffusion coefficient, tesc(p) is the mean escape time, S{p,t) represents par- 
ticle sources, and A{p) describes any additional, systematic acceleration or loss processes, 
such as shock acceleration or synchrotron/inverse-Compton emission. The quantity in square 
brackets in equation (1) describes the flux of particles through the momentum space (Tade- 
maru et al. 1971), and the source term is defined so that S{p,t) dpdt gives the number of 
particles injected into the plasma per unit volume between times t and t + dt with momenta 
between p and p + dp. The total particle number density n{t) and energy density U{t) of the 
distribution f{p,t) are computed using 

/»00 /"OO 

nit)= Anp''fip,t)dp, U{t)= Anep^fip,t)dp, (2) 

Jo Jo 

where the particle kinetic energy e is related to the Lorentz factor 7 and the particle mo- 
mentum p by 

/ 2 \ 1/2 

e = (7-l)mc2, 7=^ + 1 , (3) 
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and m and c denote the particle rest mass and the speed of hght, respectively. 

Rather than solve equation (1) directly to determine f{p,t) for a given source term 
S{p,t), it is more instructive to first solve for the Green's function, fQ{po,P,to,t), which 
satisfies the equation 



dt p^ dp \ 



Mp)fo-Dip)^ 



fa ^ ^JP - Po) K^ - to) u) 



tcsc(p) 4:Trpl 



where po is the initial momentum and to is the initial time. The source term in this equation 
corresponds to the injection of a single particle per unit volume with momentum po at time 
tg. The particle number and energy densities associated with the Green's function are given 
by 

^g(^) = / 4:Trp'^ f^{po,p,to,t)dp , U^{t) = Anep^ f^{po,p,to,t)dp . (5) 

Jo Jo 

Once the Green's function solution has been determined, the particular solution associated 

with an arbitrary source distribution S{p, t) can be computed using the integral convolution 

(e.g., Becker 2003) 

I't /"OO 

f{p,i)= / fGiPo,P,io,t)Sipo,to)dpodto , (6) 

Jo Jo 

where we have assumed that the particle injection begins at time t = and no particles are 
present in the plasma for t < 0. 



2.2. Transport CoefRcients 

In Appendix A.l, we demonstrate that for arbitrary particle energies, the mean rate of 
change of the particle momentum due to stochastic acceleration is related to the momentum 
diffusion coefficient D[p) via 



dp 
lit 



». - ?T, (^^^) ■ (^) 



The corresponding result for the mean rate of change of the kinetic energy is (see Ap- 
pendix A.2 and Miller & Ramaty 1989) 



\dt/ 



= -A (P'^D) 

stoch p^ dp 



where v is the particle speed. If the MHD wave spectrum has the power-law form W oc 
k~^ associated with Alfven and fast-mode waves, then the momentum dependences of the 
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diffusion coefficient D{p) and the mean escape time tesc{p) describing the resonant pitch-angle 
scattering of relativistic particles are given by (e.g., Dermer et al. 1996; Miller & Ramaty 
1989) 

D{p) = D, m\^ (^y , tesc(p) = U (^T^ , (9) 

\mc/ \mc/ 

where D* oc s~^ and t* oc s are constants. We shall focus on transport scenarios with g < 2, so 

that tesc is either a decreasing or constant function of the momentum p. In order to maintain 

the physical validity of the escape-timescale approach used here, we must require that tesc 

exceed the light-crossing time L/c for a source with size L. This implies a fundamental upper 

limit to the particle momentum when q < 2. 

By combining equations (7) and (9), we find that the mean rate of change of the mo- 
mentum for relativistic particles accelerated stochastically by MHD waves is given by 



dp 
'dt 



= {q + 2)D,mc(—y . (10) 

stoch \mc/ 



For simplicity, we shall assume that the momentum dependence of the additional, system- 
atic loss/acceleration processes appearing in the transport equation (1), described by the 
coefficient A{p), mimics that of the stochastic acceleration (eq. [10]). We therefore write 

A{p) = A,mc(^Y ' , (11) 



Xmc 






where the constant A^ oc s^^ determines the positive (negative) rate of systematic acceler- 
ation (losses). Note that this formulation precludes the treatment of loss processes with a 
quadratic energy dependence (e.g., inverse- Compton or synchrotron) since that would imply 
g = 3, which is outside the range considered here. However, ffist-order Fermi acceleration 
at a shock or energy losses due to Coulomb collisions can be treated by setting q = 2 with 
A* either positive or negative, respectively. This suggests that the results obtained here are 
relevant primarily for the transport of energetic ions. However, even in this application one 
needs to bear in mind that synchrotron and inverse- Compton losses will become dominant 
at sufficiently high energies (e.g., Schlickeiser 1984; Schlickeiser & Steinacker 1989; Liu et al. 
2006). 

It is convenient to transform to the dimensionless momentum and time variables x and 

y, defined by 

x = ^, y = D,t, (12) 

mc 

in which case the transport equation (4) for the Green's function becomes 

gfr , J-^ L^^, ^) - 4^ (^'^' fo)'^ ^'-' /. + "" ' "°'/'3^' '°' ■ (13) 

oy x^ ox \ ox J x^ ox ' Airm'^c^XQ 
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Xo = — , yo = D^to , (14) 

mc 

and we have introduced the dimensionless constants 

Note that x equals the particle Lorentz factor in applications involving ultrarelativistic parti- 
cles. The constant a describes the relative importance of systematic gains or losses compared 
with the stochastic process. 



2.3. Fokker-Planck Equation 

Additional physical insight can be obtained by reorganizing equation (13) in the Fokker- 
Planck form 

dN (9^ (9 

-^=g^{x''N^)- — [{q + 2 + a)x''-'N^]-9x'-'^N^ + 6{x-Xo)S{y-yo), (16) 

where we have defined the Green's function number distribution N^ using 

N^{xo,x,yo,y)=4:nm^c^x'^f^{xo,x,yo,y) ■ (17) 

The Fokker-Planck coefficients appearing in equation (16), which describe the evolution of 
the particle distribution due to the infiuence of stochastic and systematic processes, are given 
by (Reif 1965) 



(^) = (q + 2 + a)x'^-\ {U 



1 da"^ „ / dx 
^ X \ — 

2 dy ' ^dy, 

where the first coefficient describes the "broadening" of the distribution due to momentum 
space diffusion, and the second represents the mean "drift" of the particles (i.e., the average 
acceleration rate). 

Equation (16) is equivalent to equation (49) from Park & Petrosian (1995) if we set their 
parameters a^^ = 2 + a and Spp = q — 2, where Spp denotes the power-law index describing 
the momentum dependence of the escape timescale. Our particular choice for Spp refiects 
the physics of the resonant wave-particle interactions, as represented by equations (9). The 
Fokker-Planck form of equation (16) clearly reveals the fundamental nature of the transport 
process. In particular, we note that in the limit a — > 0, the drift coefficient (dx/dy) reduces to 
the purely stochastic result (eq. [10]), as expected when systematic gains/losses are excluded 



from the problem. The total number and energy densities are computed in terms of N^ 
using (cf. eq. [5]) 



where (see eq. [3]) 



N^{xo,x,yo,y)dx , U^{y) 



exN^{xo,x,yo,y)dx 



e = mc 



1/2 



{x' + iy^'-i 



(19) 



(20) 



Since the physical situation considered here corresponds to the injection of a single particle 
per unit volume at "time" y = yo, it follows that n^{yo) = 1. 



3. SOLUTION FOR THE GREEN'S FUNCTION 

The result obtained for the Green's function number distribution A''^ has two important 
special cases depending on whether q = 2 or q < 2. Park & Petrosian (1995) obtained the 
exact solution to equation (16) for the hard sphere case (Ramaty 1979) with q = 2 and 
their parameter Spp = 0, corresponding to a momentum-independent escape timescale. In 
this section we derive the exact solution to the time- dependent problem with q < 2 and 
Spp = q — 2, which describes the physics of the wave-particle interactions (see eqs. [9]). 



3.1. Laplace Transformation 

We define the Laplace transformation of A^^ using 

L{xo,x,s)= e-'y N^{xo,x,yo,y)dy . 



By operating on the Fokker-Planck equation (16) with J^ e '^^ dy, we obtain 

r/2 d 

— r (a;« L)- — \(q + 2 + a) x'^-' L\ - 9 x^-^ L - sL = -e"^^" Six - xq) , 
dx^ dx 



or, equivalent ly, 

d^L fq — 2 — a\ dL 



dx"^ 



+ 



X 



dx 



:i-g)(2 + a) 



x^ 



_e ^ 

rp2q — 2 rpQ 



L 



X'i 



(21) 



(22) 



(23) 



This equation can be transformed into standard form by introducing the new momentum 
variables 



z[x) = 



2V0 
2^q 



X 
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Zq[Xq) = 



2VO 

2^q 



X. 



2-q 



(24) 
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After some algebra, we find that the equation for L now becomes 



dPL a + 1 dL 



+ 



z^+ 



dz^ q — 2 dz 
where 



[l-q){2 + a) 
{2-qf 



sz 



4 co(2-g)2 
2v^ 



L 



CqC 



-syo 



6{z-zo) f z 



2-q 



Co 



(3-2<?)/(2-g) 



Co 



Q 



(25) 



(26) 



The solutions to equation (25) obtained for z ^ Zq that satisfy the high- and low-energy 
boundary conditions are given by 



/ (Z. Z S)- P'"/2 Ja+2)/(2-q) f ^ U{a, f3 , z) , Z > Zq 

L[Zo,z,s)-e \BM{a,P,z), z < Zo 



(27) 



where M and U denote the confluent hypergeometric functions (Abramowitz & Stegun 1970), 

and 

s + {a + 3)Ve ^ a + 3 



a = 



P 



(28) 



2(2 -g)v^ ' " 2-q 

The constants A and B appearing in equation (27) are determined by ensuring that the 
function L is continuous at z = Zq, and that it also satisfies the derivative jump condition 



, dL 
lim — — 

£-♦0 dz 



ZQ+e 



CqC 



-syo 



zo-e (2 -q)zQ \Co 



Zo 



i3^2q)/i2^q) 



-syo 



2xoV0 



(29) 



obtained by integrating the transport equation (25) with respect to z in a small range around 
the source momentum zq- 

The constant B can be eliminated by combining the continuity and derivative jump 
conditions. After some algebra, the solution obtained for A is 



A 



-syOe^o/2^(a+2)/iq-2)^^^^^^^^^ 



2xoV9W{zo) 
where W{z) denotes the Wronskian, defined by 

W{z) = M{a, P, z)^U{a, /?, z) - U{a, P, z)—M{a, p, z) . 
dz dz 



(30) 



(31) 



Using equation (13.1.22) from Abramowitz & Stegun (1970), we find that W is given by the 

exact expression 

T{p)z-^e' 



W{z) 



Via) 



(32) 
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{>Re s 



Fig. 1. — Integration contour C used to evaluate equation (37). 

which can be combined with equation (30) and the continuity relation to obtain 

r(a) ^^+(«+2)/('?-2) ^-syo g-.o/2 M{a, f3, Zo) 



A 



B 



T{a) /^+(-+^y(i'^) e-^^« e-^o/2 f/(a, p, zo) 



2xoVer{f3) 
The final solution for the Laplace transformation L can therefore be written as 



L(zo,z,s) 



where 



T{a) Zq ( z 



r(/5)2xov^ V-^o 



(a+2)/(2-q) 



(33) 
(34) 



-sy, g-(.+.o)/2 ^(^^ ^^ ^^.^) ^^^^ p^ ^^^J ^ (35) 



Zmin = ram[z, Zq) , z^^^ = max[z, Zq) . 



(36) 



3.2. Inverse Transformation 

The solution for the Green's function A^^^ can be found using the complex Mellin inversion 
integral, which states that 



1 /'7+^oo 

Na{zo,z,yo,y) = tt- e'^ L{zo,z,s) ds , 

^T^l J'y-joO 



(37) 
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where 7 is chosen so that the hne Res = 7 hes to the right of any singularities in the 
integrand. The singularities are simple poles located where |r(a;)| -^ 00, which corresponds 
to a = — n, with n = 0, 1,2, . . . Equation (28) therefore implies that there are an infinite 
number of poles situated along the real axis at s = Sn, where 

s„= [2(g-2)n-a-3] v^, n = 0,l,2,... (38) 

We can therefore employ the residue theorem to write 

/oo 
e'y L{zo, z, s) ds = 2m ^ Res(s„) , (39) 

n=0 

where C is the closed integration contour indicated in Figure 1 and Res(s„) denotes the 
residue associated with the pole at s = s„. Based on asymptotic analysis, we conclude that 
the contribution to the integral due to the curved portion of the contour vanishes in the limit 
i? — > 00, and consequently we can combine equations (37) and (39) to obtain 

00 
Ncizo,z,yo,y) = ^ Res(s„) . (40) 

ra=0 

Hence we need only evaluate the residues in order to determine the solution for the Green's 
function. 

3.3. Evaluation of the Residues 

The residues associated with the simple poles at s = s„ are evaluated using the formula 
(e.g., Butkov 1968) 

Res(s„) = hm (s — s„) e^''^ L{zq,z,s) . (41) 

Since the poles are associated with the function T{a) in equation (35), we will need to make 

use of the identity 

f— IV' 
lim {s - Sn) T{a) = ^-^ 2 (2 - g) v^ , (42) 

which follows from equations (28) and (38). Combining equations (35), (41), and (42), we 
find that the residues are given by 

Res(s„) = ^ ^ ^ .S^/^. '- - e-(^+^«)/2 Mi-n, (3, z^^) U{-n, /3, ^^^x) • 

n\V[(j)xQ \zqJ 



(43) 
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Based on equations (13.6.9) and (13.6.27) from Abramowitz & Stegun (1970), we note that 
the confluent hypergeometric functions appearing in this expression reduce to Laguerre poly- 
nomials, and therefore our result for the residue can be rewritten after some simplification 

as 



Res(s„) 



>(/3-l) 



n 



\e-My~yo)z^{2-q) /^\(-+2)/(2-.) 



T{i3 + n)xo 



zq 



-(.+.o)/2p^/3-l)(^)p^/3-l)(^^)^ (44) 



where Pn {z) denotes the Laguerre polynomial. 



3.4. Closed-Form Expression for the Green's Function 

The result for the Green's function number distribution A''^ is obtained by summing the 
residues (see eq. [40]), which yields 

(45) 
This convergent sum is a useful expression for the Green's function. However, further 
progress can be made by employing the bilinear generating function for the Laguerre poly- 
nomials, given by equation (8.976) from Gradshteyn & Ryzhik (1980). After some algebra, 
we find that the general closed-form solution can be written in the form 



A'q(xo,x, 1/0,1/) 
where 



q I X 



Xo \Xo 



(a+l)/2 






exp 



[^ + ^o)(l + ' 
2(1-0 



'/3-1 



^\/zz^i 
1-e 



Z[X) = 



2Ve 

2^q 



X 
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zo{xo) = -—— Xo 
z, q 



(3^ 



a -|- 3 
2^ 



^{y,yo) 



(46) 



Note that the solution for N^ depends on the time parameters y and yo only through the 
"age" of the injected particles, y — yo, as indicated by the form for ^. In the limit ^ ^ 0, 
corresponding to infinite escape time, the Green's function number distribution reduces to 



N^{xo,x,yo,y) 



[XXo 



i(3-g)/2 



e=o {2 - q) xl {y - yo) \xo 



X 



a/2 



— exp 



(X 
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2-q\ 



{2-qny-yo) 



'/3-1 



2(a;a;o 



l(2-g)/2 



{2-qy{y-yo) 
(48) 



The exact solution for the time-dependent Green's function describing the evolution of 
a monoenergetic initial spectrum with q < 2 given by equation (46) represents an interesting 
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new contribution to particle transport theory. The corresponding solution for the hard- 
sphere case with q = 2, given by equation (43) of Park & Petrosian (1995), can be stated in 
our notation as 



N^{xo,x,yo,y) 
where 



- ) exp 



=2 2xo^/n{y - yo) \xoJ 



— (In a; — Ina^o) 
4 (y - yo) 



(49) 



A = i^±lL + 2 + a + ^ . (50) 

We note that the general solution for A^^^ given by equation (46) agrees with equation (49) 
in the limit g -^ 2 as required. Furthermore, the equation (46) allows q to take on negative 
values if desired, and it is also applicable over a broad range of both positive and negative 
values for a. Recall that when a = 0, there are no systematic acceleration or loss pro- 
cesses included in the model. Positive (negative) values for a imply additional systematic 
acceleration (losses). 



3.5. Transition to the Stationary Solution 

The analytical solutions we have obtained for the Green's function provide a complete 
description of the response of the system to the impulsive injection of monoenergetic particles 
at any desired time. The generality of these expressions allows one to obtain the particular 
solution for the distribution function / associated with any arbitrary momentum-time source 
function S using the convolution given by equation (6). One case of special interest is the 
spectrum resulting from the continual injection of monoenergetic particles beginning at time 
t = 0, described by the source term 

S{p,t) = l^'^^P-P'^^^ ^^° ' (51) 

where Nq denotes the rate of injection of particles with momentum po per unit volume per 
unit time. We assume that no particles are present for t < 0. Combining equations (6) and 
(51), we find that the time-dependent distribution function resulting from monoenergetic 
particle injection is given by 

f{p,t)=No f f^{po,P,to,t)dto. (52) 

By transforming to the dimensionless variables x and y and employing equation (17), we 
conclude that the particular solution for the number distribution associated with continual 
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monoenergetic particle injection can be written as 

AT py 

N{x,y) =4nm^c^x^ f{x,y) = -^ N^{xo,x,yo,y) dyo , (53) 

^* Jo 

where we have used equation (14) to make the substitution dto = dyo/D^,. Since A^^ depends 
on the time parameters y and y^ only through the combination y — yo (see eqs. [46] and [49]), 
it follows that 

No '' 



N{x,y) = -f / N^{xo,x,0,y')dy' . (54) 

^* Jo 

This is a more convenient form for N{x, y) because y now appears only in the upper inte- 
gration bound. 

For general, finite values of y, the time-dependent particular solution for N{x^ y) given 
by equation (54) must be computed numerically by substituting for A^^ using either equa- 
tion (46) or (49), depending on the value of q. However, as y ^ oo, the solution rapidly 
approaches a stationary result representing a balance between injection, acceleration, and 
particle escape. The form of the stationary solution, called the steady-state Green's func- 
tion^ N^, can be obtained by directly solving the transport equation (1) with df /dt = and 
S{p, t) = Nq S{p—po)- Alternatively, the steady-state Green's function can also be computed 
by taking the limit of the time- dependent solution, which yields 

N^M = hm Nix,y) = -^ / N^ixo,x,0,y') dy' . (55) 

In the g < 2 case, we can substitute for A^^ using equation (46) and evaluate the resulting 
integral using equation (6.669.4) from Gradshteyn & Ryzhik (1980). After some algebra, we 
obtain 



2 



where P = {a + 3)/(2 — q), and 

Xmin = min(a;,a;o) , Xmax = max(a;,xo) . (57) 

Likewise, for the case with g = 2, we can substitute for A*";^ using equation (49) and then 
utilize equation (3.471.9) from Gradshteyn & Ryzhik (1980) to conclude that the steady-state 
solution is given by 

iV2(,) = -J^ (£.) '""" (^^) '^' , (58) 

5=2 2x„D.s/\\xoJ \x„,t,.J 
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where A is defined by equation (50). The steady-state solutions given by equations (56) and 
(58) agree with the results obtained by directly solving the transport equation. Due to the 
asymptotic behavior of the Bessel Ki^{z) function for large z, equation (56) indicates that 
N^ exhibits an exponential cutoff at high energies when g < 2 (Abramowitz & Stegun 1970). 
This corresponds physically to the fact that the escape timescale decreases with increasing 
particle momentum in this case. However, when q = 2 (eq. [58]), the spectrum displays a 
pure power-law behavior at high energies because the escape timescale is independent of the 
particle momentum. Specific examples illustrating these behaviors will be presented in § 4. 



3.6. Escaping Particle Distribution 

The various expressions we have obtained for the Green's function N^ describe the 
momentum distribution of the particles remaining in the plasma at time t after injection oc- 
curring at time to- However, since our model incorporates a physically realistic, momentum- 
dependent escape timescale tesc(p) given by equation (9), it is also quite interesting to com- 
pute the spectrum of the escaping particles, which may form an energetic outflow capable of 
producing observable radiation. In general, the number distribution of the escaping particles, 
N'^^'^{x,y), is related to the current distribution of particles in the plasma, N{x,y), via 

iV-^(x, y) ^ t-^l N{x, y) = OD, x'-'^ N{x, y) , (59) 

where we have used equations (9) and (15) to obtain the final result. The quantity N'^^'^{x, y) dx 
represents the number of particles escaping per unit volume per unit time with dimensionless 
momenta between x and x + dx. 

An important special case is the evolution of the escaping particle spectrum resulting 
from impulsive monoenergetic injection at dimensionless time y = y^. In this application, 
equation (59) gives the Green's function number spectrum for the escaping particles, defined 

by 

N^^\xo,x,yo,y)=eD, x'-'^ N^{xo,x,yo,y) , (60) 

where N^ is evaluated using either equation (46) or (49) depending on the value of q. We can 
use this expression for A'^^'^'^ to compute the total escaping number distribution resulting from 
an impulsive flare occurring at t = by integrating the escaping distribution with respect 
to time, which yields 

/•OO /"OO 

AAr-=(x)=/ N^^'dt = 9x'-'i N^{xo,x,0,y')dy' , (61) 

Jo Jo 

where we have used equation (60) and taken advantage of the fact that A^^ depends on y 
and I/O only through the difference y — yo (cf. eq. [54]). By comparing equations (55) and 
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(61), we deduce that 

AN-{x) = ^ x'-^ N^,{x) , (62) 

so that the total escaping spectrum is simply proportional to the steady-state Green's func- 
tion resulting from continual injection, as expected. The process considered here corresponds 
to the injection of a single particle at time t = 0, and therefore we find that the normalization 
of AN^'^ is given by 



ANl''{x)dx = l , (63) 



which provides a useful check on the numerical results. 

Another interesting example is the case of continual monoenergetic particle injection 
commencing at time t = 0, described by the source term given by equation (51). The 
time-dependent buildup of the escaping particle spectrum N'^^'^{x,y) in this scenario can be 
analyzed by using equation (54) to substitute for N{x,y) in equation (59), which yields 

N^'{x,y) = N.ex^-" r N^{xo,x,0,y')dy' . (64) 

In the limit y ^ oo, the escaping particle spectrum approaches the steady-state result 

N^^ix) = lim N^'^^ix, y) = No9 x^-'^ / NJxq, x, 0, y') dy' . (65) 

By comparing equations (55), (62), and (65), we conclude that 

N^^'{x) = OD, x^~^ N^^{x) = No AN^'^^ix) . (66) 

The final result can be combined with equation (63) to show that the escaping spectrum 
satisfies the normalization relation 

POO 

/ iV,7(a;)dx = iVo , (67) 

Jo 

as expected for the case of continual steady-state injection. 

Note that analytical expressions for the steady-state escaping spectrum N^^'^{x) can 
be obtained by substituting for N^{x) in equation (66) using either equation (56) or (58), 
depending on the value of q. We obtain 

^,;=(.) . ii^ (±) '"""^ ,.,,).-,./. j^ (fii] K._^ {^^ , ,68) 

{2-q)xQ\XQj 2 \^ 2-g y 2 y 2 - q j 

for q < 2, and 

2xoV\ \XoJ \x, ' ' 

for q = 2. 
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4. RESULTS 

The new result we have derived for the secular Green's function (eq. [46]) displays a 
rich behavior through its complex dependence on momentum, time, and the dimensionless 
parameters g, a, and 9, which are related to the fundamental physical transport coefficients 
/}*, A*, and t^, via equations (15). Here we present several example calculations in order to 
illustrate the utility of the new solution. Detailed applications to astrophysical situations, 
including active galaxies and 7-ray bursts, will be presented in subsequent papers. 

The panels on the left-hand side of Figure 2 depict the time-dependent Green's function 
solution, Nq, describing the evolution of the particle distribution in the plasma resulting 
from impulsive monoenergetic injection at y = 0. Results are presented for the hard-sphere 
scattering case (g = 2), computed using equation (49), and for the Kolmogorov (g = 5/3) 
and Kraichnan (g = 3/2) cases, evaluated using equation (46). The only acceleration mech- 
anism considered here is the stochastic acceleration associated with the second-order Fermi 
process, and therefore we set a = 0. The escape parameter 6 is set equal to unity, so that 
the timescale for escape is comparable to the diffusion timescale. As the wave turbulence 
spectrum becomes steeper (i.e., as g increases), a larger fraction of the turbulence energy is 
contained in waves with long wavelengths, which interact resonantly with higher energy par- 
ticles. Steeper turbulence spectra therefore produce harder particle distributions, as can be 
confirmed in the plots. Consequently we conclude that an ensemble of hard sphere scattering 
centers is more effective at accelerating nonthermal particles compared with a Kolmogorov 
wave spectrum, which in turn is more effective than a Kraichnan spectrum. 

The panels on the right-hand side of Figure 2 illustrate the buildup of the particle 
spectrum in the plasma, N{x,y), due to continual monoenergetic injection beginning at 
y = 0, computed by evaluating numerically the integral in equation (54). We have set a = 
and therefore the acceleration is purely stochastic. As y —^ 00, the spectrum approaches 
the steady-state form given by equation (56) for g < 2 or by equation (58) for g = 2. In 
the hard-sphere case (g = 2), the particle spectrum displays a power- law shape at high 
energies in agreement with equation (58). However, when g < 2, particle escape dominates 
over acceleration at high energies, and therefore the steady-state distribution is truncated, 
even in the absence of radiative losses. This effect produces the quasi-exponential turnovers 
exhibited by the stationary spectra when g = 5/3 and g = 3/2. Particle escape in these 
cases mimics energy losses due to, for example, synchrotron emission from electrons. 

Figures 3 and 4 illustrate the effects of varying the values of the escape parameter 6 and 
the systematic acceleration/loss parameter a for the g = 5/3 case. In Figure 3, we set a = 
and ^ = 0.1, which represents a particle escape timescale that is an order of magnitude larger 
than the corresponding case depicted in Figure 2. The longer escape time allows the particles 
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to be accelerated to higher mean energies before diffusing out of the plasma, and the decay 
of the particle density at late times takes place much more slowly than for larger values 
of 9. The hardening of the spectra causes the quasi-exponential cutoffs to move to higher 
energies, and the same effect can also be noted in the corresponding stationary solutions. 
The calculations represented in Figure 4 include additional systematic acceleration processes 
that are modeled by setting a = 2. The enhanced particle acceleration further hardens the 
spectrum and shifts the cutoff to even higher energies. Although the slope of the low-energy 
particle distribution (x < Xq) is not altered much when only 9 is varied (see Figs. 2 and 3), 
this slope becomes significantly steeper when a is increased, as indicated in Figure 4. 

In the left-hand panels of Figures 5 and 6 we plot the time-dependent solution for the 
Green's function describing the escaping particles, N'^'^, resulting from impulsive monoener- 
getic particle injection (eq. [60]) when q = 5/3. The right-hand panels illustrate the buildup 
of the escaping spectrum, N'^^'^ (eq. [64]), resulting from continual injection beginning at 
y = 0. Note the transition to the steady-state solution, N^^^ (eq. [65]), as y ^ oo. In 
Figure 5 we set ^ = 0.1 and a = 0, and in Figure 6 we set ^ = 0.1 and a = 2. Included for 
comparison are the corresponding spectra describing the particle distributions in the plasma 
at the same values of y. We point out that the escaping particle spectra are significantly 
harder than the in situ distributions, which reflects the preferential escape of the high-energy 
particles resulting from the momentum dependence of the escape timescale when q < 2 (see 
eq. [9]). 



5. DISCUSSION AND SUMMARY 

We have derived new, closed- form solutions (eqs. [46] and [68]) for the time-dependent 
Green's function representing the stochastic acceleration of relativistic ions interacting with 
MHD waves. The analytical results we have obtained describe the time-dependent distri- 
butions for both the accelerated (in situ) and the escaping particles. The Fokker-Planck 
transport equation considered here includes momentum diffusion with coefficient D{p) oc p'^, 
particle escape with mean timescale tcsc oc p'^^'^, and additional systematic acceleration/losses 
with a rate proportional to p'^"^, where p is the particle momentum and q is the index of 
the wave turbulence spectrum. This specific scenario describes the resonant interaction of 
particles with fast-mode and Alfven waves, which is one of the fundamental acceleration 
processes in high-energy astrophysics (e.g., Dermer et al. 1996). 

The new analytical result complements the work of Park & Petrosian (1995) since it is 
applicable for any q < 2 provided Spp = q — 2, where Spp is the power-law index used by these 
authors to describe the momentum dependence of the escape timescale. The most closely 
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related previous result is given by their equation (59), which treats the case with q ^ 2 
and Spp = 0, corresponding to a momentum-independent escape timescale. Our analytical 
solution (eq. [46]) agrees with theirs in the limit g — ;► 2, as expected. The general features 
of our new solution were discussed in § 4, where it was demonstrated that increasingly hard 
particle spectra result from larger values of the wave index g, smaller values of the escape 
parameter 9, and larger values of the systematic acceleration parameter a. The rich behavior 
of the solution as a function of momentum, time, and the parameters g, 9, and a provides 
useful physical insight into the nature of the coupled energetic/spatial particle transport in 
astrophysical plasmas. 

The solutions presented here can be used to describe the acceleration and transport 
of relativistic ions in astrophysical environments in which the turbulence spectrum is very 
poorly known and can be approximated by a power law, such as 7-ray bursts, active galaxies, 
magnetized coronae around black holes, and the intergalactic medium in clusters of galax- 
ies. For example, the hard X-ray emission from black-hole jet sources such as Cygnus X-1 
(Malzac et al. 2006) and the microquasar LS 5039 (Aharonian et al. 2005) could be powered 
by the stochastic acceleration of ions in a black-hole accretion disk corona that subsequently 
escape and interact with surrounding material. In this scenario, persistent acceleration of 
monoenergetic particles injected into the corona, or flaring episodes averaged over a suffi- 
ciently long time, would produce a time-averaged escaping distribution of relativistic protons 
given by equation (68) for q < 2. Assuming only stochastic acceleration, so that a = 0, the 
number distribution of escaping particles with x > xq takes the form 



iv™cc.<'-./^x,-. a^ur • ^<"^^ 



. /ii-i OC < 

2 \ 2-g ; I ^(5-2,)/2g^p/_v|^\ ^ x»to 



1/(2-9) 
l/(2-<?) 



(70) 
When the escaping hadrons collide with ambient gas or stellar wind material, they would 
generate X-rays and 7-rays via a pion production cascade with a very hard spectrum lead- 
ing up to a quasi-exponential cutoff. The Gamma ray Large Area Space Telescope (GLAST) 
telescope will be able to provide detailed spectra from Galactic black-hole sources and uniden- 
tified EGRET 7-ray sources to test for the existence of this hard component. Our new ana- 
lytical solution can also be used to model the variability of radiation from ions accelerated 
in the accretion-disk coronae of Seyfert galaxies by changing the level of turbulence. Flaring 
~ 100 MeV - GeV radiation could be weakly detected by GLAST as a consequence of this 
process. Additional applications of our work include studies of the stochastic acceleration of 
relativistic cosmic rays in 7-ray bursts (Dermer & Humi 2001). 

The results we have obtained describe both the momentum distribution of the particles 
in the plasma, and the momentum distribution of the particles that escape to form energetic 
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outflows. Since the solutions do not include inverse-Compton or synchrotron losses, which 
are usually important for energetic electrons, they are primarily applicable to cases involv- 
ing ion acceleration. In order to treat the acceleration of relativistic electrons, a generalized 
calculation including both stochastic particle acceleration and radiative losses is desirable, 
and we are currently working to extend the analytical model presented here to incorporate 
these effects. Beyond the direct utility of the new analytical solutions for probing the na- 
ture of particle acceleration in astrophysical sources, we note that they are also useful for 
benchmarking numerical simulations. 

T. L. is funded through NASA GLAST Science Investigation No. DPR-S-1563-Y. The 
work of C. D. D. and P. A. B. is supported by the Office of Naval Research. The authors 
also acknowledge the useful comments provided by the anonymous referee. 



A. APPENDIX 

In this section we explore the relationship between the momentum diffusion coefficient 
D{p) and the mean rate of change of the particle momentum and kinetic energy, using the 
integral method outlined by Subramanian et al. (1999). 



A.l. Mean Rate of Change of Momentum Due to Stochastic Acceleration 



In the case of pure stochastic acceleration, the transport equation (1) reduces to 

df _ I d 
dt p^ dp 



P D{p) - 



(Al) 



The mean momentum, (p), is defined as a function of t by 



1 f^ 
{p) = - Anp^f{p,t)dp, (A2) 

^ Jo 

where the number density n is related to / via equation (2). The associated mean rate of 
change of the momentum is given by 

Combining this relation with the transport equation (Al) yields 
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Integrating by parts once gives 



dp 
It 



n 



Anp^D—- dp 
op 



and integrating by parts a second time yields 

'|)=ir4./^(p=z))*. 

\dt / n Jq dp 



(A5) 



(A6) 



It is sufficient for our purposes to consider the evolution of the particle distribution / satis- 
fying the monoenergetic initial condition 



f{p,t) =Ao5ip-po) . 



t=to 



(A7) 



Combining this relation with equation (A6), we find that at the initial time t = to, 

1 d 



dp 
It 



t=to 



j9^ dp 



{p'D) 



(A8) 



P=PO 



Dropping the subscript "0" without loss of generality, we conclude that the mean rate of 
change of the momentum for particles with momentum p at time t is given by 



dt I p^ dp^ 



(A9) 



which agrees with equation (7). 



A. 2. Mean Rate of Change of Kinetic Energy Due to Stochastic Acceleration 

The mean kinetic energy, (e), is defined by 

(e) = -/ W^fdp, (AlO) 

where e is given by equation (3). The associated mean rate of change is given by 



/$\ = ir4,/,f,p 



\dtl 



n 



dt 



(Aii: 



which can be combined with the transport equation (Al) to obtain 



/^\ = i r Ane-(p'D^\dp 
\dt/ n Jq dp \ dp ^ 



(A12) 
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As in the previous section, we integrate by parts to find that 



\dt / n Jq dp op 






Next we employ the derivative relation 

de 



V (A14) 

dp 



and integrating by parts again to obtain 



(s)4r^^^^("^^^)*' (^1^) 



Applying the initial condition given by equation (A7), we find that at the initial time t = to 

(A16) 



/de 
\di 



= ^-^ {vp'D) 

t=to P dp^ 



p=po 



Dropping the subscript "0" without loss of generality, we conclude that the mean rate of 
change of the kinetic energy for particles with momentum p at time t is given by 

which agrees with equation (8) and with equation (13) from Miller & Ramaty (1989). 
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Fig. 2. — Green's function solutions to the time- dependent stochastic particle acceleration 
equation for xq = 1. The left panels depict the impulsive-injection solution, A''^ (eqs. [46] 
and [49]), and the right panels illustrate the response to uniform, continuous injection, N 
(eq. [54]), with Nq = D^, = 1. Note that N approaches the corresponding steady-state 
solution (eqs. [56] and [58]) as y increases. The indices of the wave turbulence spectra are 
indicated, with q = 2 for hard-sphere scattering, q = 5/3 for a Kolmogorov cascade, and 
q = 3/2 for a Kraichnan cascade. 
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Fig. 3. — Evolution of the particle distribution in the plasma resulting from monoenergetic 
injection with q = 5/3, a = 0, 9 = 0.1, and Xq = 1. Panel (a) depicts the Green's function, 
N^, resulting from impulsive monoenergetic injection (eq. [46]), and panel (b) illustrates the 
response to continual monoenergetic injection (eq. [54]) and the corresponding steady-state 
solution (eq. [56]) for A^o = -D* = 1. 
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Fig. 4. — Same as Fig. 3, except a = 2. 
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Fig. 5. — Evolution of the particle distribution resulting from monoenergetic injection with 
q = 5/3, a = 0, 9 = 0.1, xq = 1, Nq = 1, and D* = 1. Panel (a) treats the case of impulsive 
injection, with the thin lines representing the particle distribution in the plasma (eq. [46]) 
and the heavy lines denoting the escaping particle spectrum (eq. [60]) in units of 6. Panel (6) 
illustrates the response to continual monoenergetic injection for the particle distribution in 
the plasma (eqs. [54] and [56]; thin lines) and for the escaping particle spectrum (eqs. [64] 
and [66]; heavy lines). See the discussion in the text. 
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Fig. 6. — Same as Fig. 5, except a = 2. 



